home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / External Modules / External modules sources / C / Contour plotting / ContourPlotter.c < prev    next >
C/C++ Source or Header  |  1996-04-21  |  12KB  |  410 lines
  1. #include "proFit_interface.h"
  2.  
  3. #ifndef __CONTOUR_PLOTTER__
  4. #include "ContourPlotter.h"
  5. #endif
  6.  
  7.  
  8. typedef struct MyPrivatePlottingData
  9. {
  10.     long            nrRows, nrCols;            // the number of rows and columns to scan
  11.     void**            mesh;                    // pointer to a 2-d array of mesh points, see below
  12.     Boolean            exit;                    // set to true if we must exit
  13.     double            cMin, cMax;                // values corresponding to first and last column
  14.     double            rMin, rMax;                // values corresponding to first and last row
  15.     double            rUnit, cUnit;            // (xiMax-xiMin)/(nrRows(nrCols)-1)
  16.  
  17. }MyPlottingData;
  18.  
  19.  
  20. // the following struct gives info about a single mesh point
  21. typedef struct
  22. {
  23.     double    value;                    /* the value at the given point */
  24.     Boolean    valid;                    /* true if this point has valid data */
  25.     Boolean    atRight;                /* true if a contour passes the mesh to the right of the point */
  26.     Boolean    atBottom;                /* true if a contour passes the mess to the bottom of the point */
  27. } MeshPoint;
  28.  
  29.  
  30.  
  31.  
  32. /* the following is the size of an entry in the mesh array - note that we make sure */
  33. /* that the size is a multiple of 8 for increasing speed on systems where doubles should not cross */
  34. /* 2- or 4 byte boundaries */
  35. #define meshPointSize ((long)((sizeof(MeshPoint)+7)&0xFFFFFF8))
  36.  
  37. #define ABS(x)  ((x)<0? (-(x)): (x))
  38.  
  39.  
  40. static MeshPoint* GetMeshPoint (long row, long column, MyPlottingData* const v)
  41.     // returns the address of the mesh point at the given row/column
  42.     // row: 0..v->nrRows-1, col: 0..v->nrCols-1
  43. {
  44.     long    offset = meshPointSize * (row*v->nrCols + column);
  45.  
  46.     return (MeshPoint*)(*(char**)v->mesh + offset);
  47. }
  48.  
  49.  
  50. static Boolean TestIfExit(MyPlottingData* const v)
  51.     /* returns true if v->exit is true or if user pressed cmd-. */
  52. {
  53.     if (TestStop()) v->exit = true;
  54.     return v->exit;
  55. }
  56.  
  57. static void AddToArray(double** array, double value, long nrValues)
  58.     /* add value to array, array presently holding nrValues values */
  59.     /* resizging the array if necessary */
  60. {
  61.     long    requiredSize = (nrValues+1)*sizeof(double);
  62.     if (GetHandleSize((Handle)array) < requiredSize)
  63.     {    SetHandleSize((Handle)array, requiredSize+1024);
  64.         if (MemError() != noErr) return;
  65.     }
  66.     (*array)[nrValues] = value;
  67. }
  68.  
  69. static void Scale(double row, double col, double* const x1, double* const x2, MyPlottingData* const v)
  70.     // converts units in rows/columns to units in x1, x2
  71. {
  72.     *x1 = v->cMin + col * v->cUnit;
  73.     *x2 = v->rMin + row * v->rUnit;
  74. }
  75.  
  76. static void DrawOneContour(double level, long startR, long startC, MyPlottingData* const v)
  77.     /* draw one contour starting at row/col startR/startC */
  78. {
  79.     long        r, c;
  80.     Boolean        done = false;
  81.     int            where;                                /* 0 if present point is at right, 1 if at button */
  82.     MeshPoint*    info = GetMeshPoint(startR, startC, v);        /* the point that we are currently looking at */
  83.     MeshPoint*    origInfo = info;                            /* the point where we started */
  84.     double        **xValues1=nil, **yValues1=nil;        /* x- and y-values at first branch */
  85.     double        **xValues2=nil, **yValues2=nil;        /* x- and y-values at second branch */
  86.     double        **xValues, **yValues;
  87.     long        nrValues1=0, nrValues2=0;            /* number of values in first and second branch */
  88.     long        *nrValues;
  89.     int            branch=0;                            /* 0: running first branch, 1: second branch */
  90.     int            firstWhere;                            /* the value of "where" at the first point */
  91.  
  92.     xValues1 = (double**)NewHandle(256);
  93.     xValues2 = (double**)NewHandle(256);
  94.     yValues1 = (double**)NewHandle(256);
  95.     yValues2 = (double**)NewHandle(256);
  96.     if (xValues1==nil || xValues2==nil || yValues1==nil || yValues2==nil) goto done;
  97.  
  98.  
  99.     firstWhere = !info->atRight;
  100.     while (!TestIfExit(v) && branch < 2)            /* while there's something to trace */
  101.     {
  102.         if (branch == 0)
  103.         {    where = firstWhere;
  104.             xValues = xValues1;
  105.             yValues = yValues1;
  106.             nrValues = &nrValues1;
  107.         }
  108.         else                                        /* if second branch */
  109.         {    where = firstWhere;
  110.             xValues = xValues2;
  111.             yValues = yValues2;
  112.             nrValues = &nrValues2;
  113.         }
  114.         r=startR, c=startC;
  115.         info = GetMeshPoint(r, c, v);
  116.         branch += 1;
  117.         while (!TestIfExit(v))                        /* trace one branch */
  118.         {
  119.             MeshPoint*    next;
  120.             double        d;
  121.     
  122.             if (branch == 1 && *nrValues == 2)        /* after two steps */
  123.             {    if (firstWhere == 0)                /* we must set original point back to true (if loop is closed, original point is used twice) */
  124.                     origInfo->atRight = true;
  125.                 else
  126.                     origInfo->atBottom = true;
  127.             }
  128.  
  129.             if (where==0)                            /* if we are to the right of the present mesh point */
  130.             {
  131.                 info->atRight = false;                /* delete this point - it is traced now */
  132.                 next = GetMeshPoint(r, c+1, v);
  133.                 if (next->valid==false)
  134.                     break;
  135.     
  136.                 d = ABS(next->value - info->value);
  137.                 if (d == 0)
  138.                     break;
  139.                 
  140.                 AddToArray(xValues, r, *nrValues);
  141.                 AddToArray(yValues, c + ABS(info->value-level)/d, (*nrValues)++);
  142.     
  143.                 if (info->atBottom)
  144.                 {    where = 1;
  145.                     continue;
  146.                 }
  147.  
  148.                 if (c+1 != v->nrCols)                    /* if not yet at right of grid */
  149.                 {    next = GetMeshPoint(r, c+1, v);
  150.                     if (next->valid && next->atBottom)
  151.                     {    c += 1;
  152.                         where = 1;
  153.                         info = next;
  154.                         continue;
  155.                     }
  156.                 }
  157.                 if (r+1 != v->nrRows)                    /* if not at bottom of grid */
  158.                 {    next = GetMeshPoint(r+1, c, v);
  159.                     if (next->valid && next->atRight)
  160.                     {    r += 1;
  161.                         where = 0;
  162.                         info = next;
  163.                         continue;
  164.                     }
  165.                 }
  166.                 if (r != 0)                                /* if not at top of grid */
  167.                 {
  168.                     next = GetMeshPoint(r-1, c, v);
  169.                     if (next->valid && (next->atRight || next->atBottom))
  170.                     {    info = next;
  171.                         r -= 1;
  172.                         where = !info->atRight;
  173.                         continue;
  174.                     }
  175.                     if (c+1 != v->nrCols)                /* if not yet at right of grid */
  176.                     {    next = GetMeshPoint(r-1, c+1, v);
  177.                         if (next->valid && next->atBottom)
  178.                         {    r -= 1;
  179.                             c += 1;
  180.                             where = 1;
  181.                             info = next;
  182.                             continue;
  183.                         }
  184.                     }
  185.                 }
  186.                 break;                                    /* did not find continuation */
  187.             }//if where==0
  188.             else
  189.             {
  190.                 info->atBottom = false;                    /* delete this point - it is traced now */
  191.                 next = GetMeshPoint(r+1, c, v);
  192.                 if (next->valid==false)
  193.                     break;
  194.     
  195.                 d = ABS(next->value - info->value);
  196.                 if (d == 0)
  197.                     break;
  198.                 
  199.                 AddToArray(xValues, r + ABS(info->value-level)/d, (*nrValues));
  200.                 AddToArray(yValues, c, (*nrValues)++);
  201.     
  202.                 if (info->atRight)
  203.                 {    where = 0;
  204.                     continue;
  205.                 }
  206.  
  207.                 if (r+1 != v->nrRows)                    /* if not yet at bottom of grid */
  208.                 {    next = GetMeshPoint(r+1, c, v);
  209.                     if (next->valid && next->atRight)
  210.                     {        r += 1;
  211.                             where = 0;
  212.                             info = next;
  213.                             continue;
  214.                     }
  215.                 }
  216.                 if (c+1 != v->nrCols)                    /* if not yet at right of grid */
  217.                 {    next = GetMeshPoint(r, c+1, v);
  218.                     if (next->valid && next->atBottom)
  219.                     {    c += 1;
  220.                         where = 1;
  221.                         info = next;
  222.                         continue;
  223.                     }
  224.                 }
  225.                 if (c != 0)                                // if not at left of grid
  226.                 {    next = GetMeshPoint(r, c-1, v);
  227.                     if (next->valid && (next->atRight || next->atBottom))
  228.                     {    info = next;
  229.                         c -= 1;
  230.                         where = !info->atRight;
  231.                         continue;
  232.                     }
  233.                     if (r+1 != v->nrRows)                /* if not yet at bottom of grid */
  234.                     {    next = GetMeshPoint(r+1, c-1, v);
  235.                         if (next->valid && next->atRight)
  236.                         {    c -= 1;
  237.                             r += 1;
  238.                             where = 0;
  239.                             info = next;
  240.                             continue;
  241.                         }
  242.                     }
  243.                 }
  244.                 break;                                    /* did not find continuation */
  245.             }//if where==0
  246.         }//while
  247.     }//while
  248.     if (TestIfExit(v)) goto done;
  249.     {
  250.         Boolean    first=true;
  251.         long    i;
  252.         double    x, y, xLast, yLast;
  253.  
  254.         if (nrValues1>0)
  255.         {    for (i=nrValues1-1; i>=0; i--)
  256.             {    Scale((*xValues1)[i], (*yValues1)[i], &x, &y, v);
  257.                 if (first) GrMoveTo(x, y);
  258.                 else GrLineTo(x, y);
  259.                 first = false;
  260.             }
  261.             xLast = x;
  262.             yLast = y;
  263.         }
  264.         for (i=0; i < nrValues2; i++)
  265.         {    Scale((*xValues2)[i], (*yValues2)[i], &x, &y, v);
  266.             if (nrValues1 > 0 && i == 0)                /* if the first point of this branch */
  267.             {    if (x==xLast && y==yLast)                /* if last branch ended with same point as */
  268.                     continue;                            /* this point starts with */
  269.             }
  270.             if (first) GrMoveTo(x, y);
  271.             else GrLineTo(x, y);
  272.             first = false;
  273.         }
  274.     }
  275.  
  276. done:
  277.     if (xValues1 != nil) DisposeHandle((Handle)xValues1);
  278.     if (xValues2 != nil) DisposeHandle((Handle)xValues2);
  279.     if (yValues1 != nil) DisposeHandle((Handle)yValues1);
  280.     if (yValues2 != nil) DisposeHandle((Handle)yValues2);
  281. }
  282.  
  283.  
  284.  
  285. PlottingData* InitContourPlotter(long nrRows, long nrCols,
  286.         double x1Min, double x1Max, double x2Min, double x2Max)
  287.     // call this routine to initialize the contour plotter for 
  288.     // working on a matrix of size nrRows, nrCols
  289.     // x1Min, x1Max, etc. correspond to the numerical values of the bounds of the matrix
  290.     // returns nil if not enough memory
  291.     // if InitContour is successful, it returns the above defined struct.
  292. {
  293.  
  294.     MyPlottingData* v = (MyPlottingData*)NewPtrClear(sizeof(MyPlottingData));
  295.     long    size;
  296.     OSErr    err;
  297.     Handle    myMemory;
  298.  
  299.     if (v==nil) return nil;
  300.  
  301.     v->exit = false;
  302.     v->mesh = nil;
  303.     v->nrRows = nrRows;
  304.     v->nrCols = nrCols;
  305.     v->cMin = x1Min; v->cMax = x1Max;
  306.     v->rMin = x2Min; v->rMax = x2Max;
  307.  
  308.     v->rUnit = (x2Max - x2Min)/(nrRows-1);
  309.     v->cUnit = (x1Max - x1Min)/(nrCols-1);
  310.  
  311.     size = nrRows*nrCols*meshPointSize;
  312.     
  313.     if (TempFreeMem() > size+50000)                            // if there's a lot of temporary mem
  314.         myMemory = TempNewHandle(size, &err);                // try using temporary memory
  315.     if (myMemory == nil)
  316.         myMemory = NewHandle(size);                            // try using heap memory
  317.     if (myMemory == nil)                                    // if allocation failed
  318.     {    DisposePtr((Ptr)v);
  319.         return nil;
  320.     }
  321.     HLock(myMemory);
  322.     v->mesh = (void**)myMemory;
  323.     return (PlottingData*)v;
  324. }//InitContourPlotter
  325.  
  326. void DisposeContourPlotter(PlottingData* p)
  327.     // deallocates all memory allocated by InitContourPlotter    
  328.     // does nothing if v==nil
  329. {
  330.     MyPlottingData* v = (MyPlottingData*)p;
  331.     if (v)
  332.     {    if (v->mesh) DisposeHandle((Handle)v->mesh);
  333.         DisposePtr((Ptr)v);
  334.     }
  335. }
  336.  
  337.  
  338. Boolean CalculateContourMatrix(FunctionPtr func, void* param, PlottingData* p)
  339.     // after having called InitContourPlotter, call this routine to initialize
  340.     // the plotting matrix.
  341.     // func: a function that returns x3 for each matrix point x1, x2
  342.     // param is passed to func
  343. {
  344.     double    x1, x2;
  345.     long    r, c;
  346.     MyPlottingData* v = (MyPlottingData*)p;
  347.  
  348.     x2 = v->rMin;
  349.     for (r = 0; r < v->nrRows; r++)            // initialize the mesh point array
  350.     {
  351.         x1 = v->cMin;
  352.         for (c = 0; c < v->nrCols; c++)
  353.         {    MeshPoint*    mesh = GetMeshPoint(r, c, v);
  354.             mesh->valid = (*func)(x1, x2, &mesh->value, param);
  355.             if (TestIfExit(v)) return false;
  356.             x1 += v->cUnit;
  357.         }
  358.         x2 += v->rUnit;
  359.     }
  360.     return true;
  361. }
  362.  
  363.  
  364.  
  365. Boolean DrawContour(double level, PlottingData* p)
  366.     // after having called CalculateContourMatrix, call this routine for all levels
  367.     // (i.e. x3-values) where you want to have a contour line
  368.     // Returns false if execution must be stopped
  369.  
  370. {
  371.     long r, c;
  372.     MyPlottingData* v = (MyPlottingData*)p;
  373.  
  374.     // first check where we have crossing points in the matrix:
  375.     for (r = 0; r < v->nrRows && !TestIfExit(v); r++)
  376.     {
  377.         for (c = 0; c < v->nrCols && !TestIfExit(v); c++)
  378.         {
  379.             MeshPoint*    mesh = GetMeshPoint(r, c, v);
  380.             if (mesh->valid)                                    // if there's valid data in this cell
  381.             {
  382.                 Boolean    thisIsBelowLevel;
  383.                 mesh->atRight = mesh->atBottom = false;
  384.                 thisIsBelowLevel = mesh->value < level;
  385.                 if (c < v->nrCols-1)
  386.                 {    MeshPoint* right = GetMeshPoint(r, c+1, v);        // this is not very efficient, might use buffering here
  387.                     if (right->valid && (thisIsBelowLevel != (right->value < level)))
  388.                         mesh->atRight = true;
  389.                 }
  390.                 if (r < v->nrRows-1)
  391.                 {    MeshPoint* bottom = GetMeshPoint(r+1, c, v);
  392.                     if (bottom->valid && (thisIsBelowLevel != (bottom->value < level)))
  393.                         mesh->atBottom = true;
  394.                 }
  395.             }
  396.         }
  397.     }//for r
  398.  
  399.     // now follow all crossing points to create the contours:
  400.     for (r = 0; r < v->nrRows-1 && !TestIfExit(v); r++)
  401.     {    for (c = 0; c < v->nrCols-1 && !TestIfExit(v); c++)
  402.         {
  403.             MeshPoint*    info = GetMeshPoint(r, c, v);
  404.             if (info->atRight) DrawOneContour(level, r, c, v);
  405.             if (info->atBottom) DrawOneContour(level, r, c, v);
  406.         } /*for c*/
  407.     } /*for r*/
  408.     return !TestIfExit(v);
  409. }
  410.